ca_dgw <- read_sf(here("ca_dgw"),
layer = "F2013_DBGS_Points_20150720_093252") %>%
clean_names()
st_crs(ca_dgw)
## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]
ca_counties <- read_sf(here("ca_counties"),
layer = "CA_Counties_TIGER2016") %>%
clean_names() %>%
select(name)
st_crs(ca_counties)
## Coordinate Reference System:
## User input: WGS 84 / Pseudo-Mercator
## wkt:
## PROJCRS["WGS 84 / Pseudo-Mercator",
## BASEGEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]],
## CONVERSION["Popular Visualisation Pseudo-Mercator",
## METHOD["Popular Visualisation Pseudo Mercator",
## ID["EPSG",1024]],
## PARAMETER["Latitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8801]],
## PARAMETER["Longitude of natural origin",0,
## ANGLEUNIT["degree",0.0174532925199433],
## ID["EPSG",8802]],
## PARAMETER["False easting",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8806]],
## PARAMETER["False northing",0,
## LENGTHUNIT["metre",1],
## ID["EPSG",8807]]],
## CS[Cartesian,2],
## AXIS["easting (X)",east,
## ORDER[1],
## LENGTHUNIT["metre",1]],
## AXIS["northing (Y)",north,
## ORDER[2],
## LENGTHUNIT["metre",1]],
## USAGE[
## SCOPE["unknown"],
## AREA["World - 85°S to 85°N"],
## BBOX[-85.06,-180,85.06,180]],
## ID["EPSG",3857]]
ca_counties <- st_transform(ca_counties, st_crs(ca_dgw))
st_crs(ca_counties)
## Coordinate Reference System:
## User input: WGS 84
## wkt:
## GEOGCRS["WGS 84",
## DATUM["World Geodetic System 1984",
## ELLIPSOID["WGS 84",6378137,298.257223563,
## LENGTHUNIT["metre",1]]],
## PRIMEM["Greenwich",0,
## ANGLEUNIT["degree",0.0174532925199433]],
## CS[ellipsoidal,2],
## AXIS["latitude",north,
## ORDER[1],
## ANGLEUNIT["degree",0.0174532925199433]],
## AXIS["longitude",east,
## ORDER[2],
## ANGLEUNIT["degree",0.0174532925199433]],
## ID["EPSG",4326]]
ggplot() +
geom_sf(data = ca_counties) +
geom_sf(data = ca_dgw, aes(color = dgbs))
tmap_mode("view")
## tmap mode set to interactive viewing
tm_shape(ca_dgw) +
tm_dots("dgbs")
## Variable(s) "dgbs" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.
sj_county <- ca_counties %>%
filter(name == "San Joaquin")
sj_depth <- ca_dgw %>%
st_intersection(sj_county)
## although coordinates are longitude/latitude, st_intersection assumes that they are planar
## Warning: attribute variables are assumed to be spatially constant throughout all
## geometries
plot(sj_depth)
## Warning: plotting the first 10 out of 18 attributes; use max.plot = 18 to plot
## all
plot(sj_county)
ggplot() +
geom_sf(data = sj_county) +
geom_sf(data = sj_depth, aes(color = dgbs))
# Find dupes
well_duplicates <- sj_depth %>%
get_dupes(latitude, longitude)
# Get rid of dupes
sj_depth <- sj_depth %>%
filter(!local_well %in% well_duplicates$local_well)
sj_depth %>%
get_dupes(latitude, longitude)
## No duplicate combinations found of: latitude, longitude
## Simple feature collection with 0 features and 19 fields
## bbox: xmin: NA ymin: NA xmax: NA ymax: NA
## geographic CRS: WGS 84
## # A tibble: 0 x 20
## # … with 20 variables: latitude <dbl>, longitude <dbl>, dupe_count <int>,
## # site_code <chr>, local_well <chr>, state_well <chr>, wcr_number <chr>,
## # well_use <dbl>, msmt_date <date>, msmt_agenc <dbl>, wsel <dbl>, dgbs <dbl>,
## # rp_elevati <dbl>, gs_elevati <dbl>, msmt_metho <dbl>, msmt_issue <dbl>,
## # msmt_comme <chr>, link_to_wd <chr>, name <chr>, geometry <GEOMETRY [°]>
Create a variogram
sj_dgw_vgm <- variogram(dgbs ~ 1, data = sj_depth)
plot(sj_dgw_vgm)
# Find model that fits spatial variogram
sj_dgw_vgm_fit <- fit.variogram(sj_dgw_vgm,
model = vgm(nugget = 20,
psill = 3000,
range = 30,
model = "Gau"))
sj_dgw_vgm_fit
## model psill range
## 1 Nug 102.3052 0.00000
## 2 Gau 2843.6996 17.18188
plot(sj_dgw_vgm, sj_dgw_vgm_fit)
Spatial kriging (interpolation)
sj_grid <- st_bbox(sj_county) %>%
st_as_stars(dx = 0.01, dy = 0.01) %>%
st_set_crs(4326) %>%
st_crop(sj_county)
## although coordinates are longitude/latitude, st_union assumes that they are planar
## although coordinates are longitude/latitude, st_intersects assumes that they are planar
plot(sj_grid)
sj_dgw_krige <- krige(dgbs ~ 1, sj_depth, sj_grid,
model = sj_dgw_vgm_fit)
## Warning in proj4string(d$data): CRS object has comment, which is lost in output
## Warning in proj4string(newdata): CRS object has comment, which is lost in output
## [using ordinary kriging]
## Warning in proj4string(newdata): CRS object has comment, which is lost in output
## Warning in sp::proj4string(.x): CRS object has comment, which is lost in output
## Warning in sp::proj4string(.x): CRS object has comment, which is lost in output
plot(sj_dgw_krige)
